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Abstract 

A  combined  thermal  and  simple  kinetic  model  is  applied  to  a  small  single  solid  oxide  fuel  cell  (SOFC,  20  cm2  square  cell  of  anode 
supported  electrolyte  with  1  cm2  active  area).  The  aim  was  to  compute  the  temperature  field  to  analyze  the  operating  conditions  of  the  active 
area  of  the  cell.  The  cell  is  simulated  in  the  conditions  used  for  electrochemical  characterization,  i.e.  at  negligible  fuel  utilization  and 
thermally  non- adiabatic. 

Though  using  a  very  simple  kinetic  model,  the  simulated  polarization  curves  fit  the  experimental  results  at  high  current  density  where  the 
curve  slope  tends  to  decrease.  The  simulated  temperature  field  shows  clearly  that  this  aspect  of  the  curve  is  explained  by  a  local  temperature 
increase  (of  30  K  at  1  A/cm2).  The  temperature  of  the  active  area  can  be  plotted  versus  the  current  and  this  result  fits  previous  measurements. 

Finally,  the  temperature  profile  simulated  shows  that,  with  the  usual  temperature  measurement  used  in  standard  electrochemical  testing,  a 
few  millimeters  away  from  the  active  surface,  does  not  detect  the  effective  cell  temperature.  In  a  stack  modeling  perspective,  the  quality  of  the 
kinetic  model  used  is  essential.  This  model,  combined  with  a  parameter  estimation  tool,  allows  using  experimental  results  to  end,  from  routine 
measurements,  with  an  accurate  and  up-to-date  kinetic  model. 

©  2003  Elsevier  Science  B.V.  All  rights  reserved. 

Keywords:  Solid  oxide  fuel  cell;  Model;  Kinetics;  Local  heating 


1.  Introduction 

In  operation,  the  temperature  field  in  a  solid  oxide  fuel  cell 
(SOFC)  stack  can  show  gradients  up  to  100  K  [1-3].  As  the 
electrochemical  reactions  are  activated  with  temperature,  a 
non-homogeneous  current  density  field  is  caused  by  fuel  and 
air  utilization  and  local  temperature  increase.  Kinetics  of  the 
electrochemical  reaction  need  a  special  focus  as  they  are 
responsible  both  for  the  power  output  and  the  heat  sources  of 
the  stack.  Reports  are  published  on  SOFC  electrode  kinetics 
treating  in  detail  the  reaction  path  [4,5],  micro  structure  of 
the  anode  cermet  [6]  and  cathode.  For  stack  modeling,  this 
information  is  too  complex,  and  mostly  applies  only  at  low 
overpotential  and  current.  Furthermore,  the  results  published 
are  valid  for  the  particular  materials  tested,  which  usually 
differ  from  the  ones  prepared  in-house  for  electrochemical 
testing  and  repeat  element  assembling  [7]. 

This  work  intends  to  present  a  way  of  using  routine 
electrochemical  results  to  determine  parameters  for  a  kinetic 
model.  This  could  allow  to  update  a  kinetic  model  following 
continuous  cell  development.  The  usual  electrochemical 
charaterization  is  simple:  a  single  cell  (anode  supported 
electrolyte)  with  1  cm2  cathode  is  tested  at  750  or  800  °C 
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(oven  temperature).  For  the  test,  a  thermocouple  is  usually 
placed  on  the  cathode  side  a  few  millimeters  away  from  the 
active  area;  this  temperature  rises  a  little  with  current 
(around  4  K  at  1.5  A/cm2).  However,  there  is  evidence  that 
the  temperature  increase  on  the  active  surface  can  reach 
considerably  higher  values  (around  20-25  °C  at  1  A/cm2 
reported  [8]). 

To  use  such  routine  electrochemical  testing  for  a  repeat 
element  model,  a  specific  model  representing  the  thermal 
phenomenon  on  a  single  cell  has  been  developed.  The 
temperature  versus  current  and  the  current-potential  char¬ 
acteristics  could  be  reproduced. 

2.  Methodology 

This  section  presents  the  method  used  for  calibration  of 
the  model  developed  with  the  experiments. 

2.7.  Model  and  experiment 

The  experiment  carried  out  differs  from  the  routine 
electrochemical  test.  Two  thermocouples  were  employed: 
one  on  the  usual  site  (a  few  millimeters  away  from  the 
cathode)  and  one  on  the  cathode  itself  to  monitor  the  local 
temperature  of  the  active  area. 
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The  temperature  increase  on  the  active  surface  of  the  cell 
is  obviously  caused  by  the  electrochemical  processes  tak¬ 
ing  place  on  the  electrodes.  However,  the  fact  that  tem¬ 
perature  measurement  close  to  that  surface  is  unable  to 
detect  this  phenomenon  is  explained  by  the  thermal  proper¬ 
ties  and  thermal  boundary  conditions  of  the  cell.  A  thermal 
model  describing  the  reactions,  heat  transfer  and  energy 
balance  phenomena  has  been  developed,  and  is  presented  in 
Section  3. 


2.2.  Parameter  estimation:  problem  formulation 


Model  and  experiments  are  compared  using  a  parameter 
estimation  method.  This  methodology  is  well  known  in 
process  modeling  and  is  applied  here  to  determine  para¬ 
meters  for  the  thermal  model  and  for  the  kinetic  model.  The 
general  formulation  for  this  kind  of  problem  is  described 
here.  A  model  can  take  the  form  of  a  set  of  equations  which 
must  be  satisfied: 


F(X,zA0)  =  0 


where  X  are  the  non-measured  variables  and  z  the  measured 
variables.  The  set  of  parameters  6  is  divided  in  parameters  8 
that  will  be  identified  by  the  experiment  and  6  which  are 
fixed.  The  algorithm  used  is  the  tool  for  parameter  estima¬ 
tion  in  the  gPROMS  [9,10]  package  and  the  objective 
function  f  to  be  minimized  is  the  following: 
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where  N  is  the  number  of  measurements;  9  the  set  of 
parameters  to  be  optimised;  NE  the  number  of  experiments 
performed;  NV;  the  number  of  variables  in  the  ith  experi¬ 
ment;  NM ij  the  number  of  measurements  of  the  jth  variable 
in  the  ith  experiment;  Zijk  is  the  kth  measurement  value  of 
variable  j  in  experiment  i;  Zijk  is  variance  of  the  kth  measure¬ 
ment  of  variable  j  in  experiment  i;  and  Zijk  the  kth  model 
predicted  value  of  variable  j  in  experiment  i. 

Experiments  with  temperature,  current  and  potential  mea¬ 
surements  are  used  to  calibrate  the  thermal  model  (i.e.  the 
temperature  response  to  the  current  variation).  Temperature 
is  measured  on  the  active  surface  and  a  few  millimeters  away 
from  it.  The  second  step  is  to  verify  that  with  a  routine 
electrochemical  experiment,  values  for  the  kinetics  can  be 
found  and  used  with  enough  confidence. 

The  difference  between  measured  and  predicted  values  is 
minimized  by  the  algorithm.  Optimization  is  required  due  to 
the  non-linearity  and  complexity  of  the  model  which  is 
going  to  be  presented  further. 

The  methodology  is  based  on  a  two-step  optimization. 
First,  experiments  with  several  temperature  measurements 
are  made  to  validate  the  model  and  determine  the  values  for 
some  parameters  into  the  thermal  model  itself.  Once  those 


parameters  are  identified,  the  tool  is  ready  to  be  used  on  the 
routine  electrochemical  characterization  to  determine  para¬ 
meters  for  the  kinetic  model. 


3.  Thermal  model 

3.1.  Experimental  set-up  and  model  description 

The  thermal  model  intends  to  represent  a  cell  of  small  size 
(16-25  cm2)  with  a  1  cm2  screen-printed  cathode  (i.e.  active 
area),  placed  into  a  seal-less  set-up  consisting  of  two  spring- 
loaded  flanges  with  a  single  gas  inlet  each.  Current  collec¬ 
tion  is  achieved  by  pressing  nickel  mesh  on  the  anode  and 
platinum  grid  on  the  cathode.  The  cell  is  kept  between  two 
zirconia  felts  (90%  porous).  The  oven  is  set  to  temperature  in 
the  range  between  650  and  850  °C.  Current-potential  mea¬ 
surements  were  recorded  with  a  potentiostat.  Two  thin 
thermocouples  of  0. 1  mm  diameter  (type  K)  have  been  used 
for  the  temperature  measurements  (see  Figs.  1  and  2). 

As  the  cell  has  air  and  fuel  flowing  on  each  side,  con¬ 
vective  heat  transfer  takes  place.  Considering  the  flanges, 
heat  exchange  is  mainly  due  to  radiation  exchange  with  the 
oven’s  walls.  The  thermal  model  represents  the  following 
elements: 


radiative  loss 


Fig.  1.  Seal-less  set-up  using  flanges  and  zirconia  oxide  felt  to  hold  the 
cell  in  place  and  feed  the  gases. 
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screen  printed 
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Fig.  2.  Scheme  of  the  cell  modeled. 
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•  the  flanges  on  air  and  fuel  side; 

•  the  cell  (as  a  single  solid  layer); 

•  the  gases  between  the  cell  and  the  flanges. 

Each  of  these  elements  is  described  with  an  energy 
balance  equation  computing  the  temperature  field.  The  fuel 
and  oxygen  utilization  in  these  tests  is  inferior  to  10%  and 
therefore  mass  balance  equations  for  the  gases  are  not 
computed.  Introduction  of  this  set  of  mass  balance  equations 
leads  to  convergence  problems  because  of  the  discontinuity 
between  the  active  area  and  the  rest  of  the  cell.  Finally  the 
assumptions  leading  to  the  energy  balance  equations  are: 

•  thermal  boundary  conditions  on  the  flanges  are  non- 
adiabatic,  the  losses  are  essentially  radiative  (horizontal 
mounting  limits  natural  convection); 

•  a  simplified  description  of  the  flow  pattern  in  2D  is  made 
with  potential  flow  [11]; 

•  post-combustion  of  the  unreacted  hydrogen  is  assumed  to 
take  place  at  the  edge  of  the  cell  to  the  extent  of  available 
excess  air; 

•  convective  heat  transfer  is  dominating  the  exchange 
between  solids  and  gases,  the  flow  is  considered  laminar 
and  fully  developed  as  first  approximation; 

•  kinetic  description  is  kept  simple  (ohmic  losses  and 
activation  term)  to  be  easily  handled  in  a  future  stack 
model; 

•  the  temperature  gradient  in  the  ^-direction  of  the  cell  is 
neglected. 

The  system  boundary  of  the  model  is  indicated  in  Fig.  1 
and  include  the  end  of  the  gas  feed  tubes. 

The  thermal  properties  for  anode  supported  cells  are  quite 
different  from  those  of  the  electrolyte  supported  cell  where  the 
yttria-stabilized  zirconia  (YSZ)  dominates  the  heat  conduc¬ 
tion  [12].  The  thermal  properties  of  the  cermet  Ni/YSZ  are 
very  sensitive  to  the  composition  and  porosity.  Between  32 
and  44  vol.%  Ni,  and  the  corresponding  porosities,  the  cell 
conductivity  varies  between  4.4  and  10  W/m  K.  These  values 
are  higher  than  the  Y SZ  thermal  conductivity  which  is  usually 
reported  to  be  on  the  order  of  2  W/(m  K). 

3.2.  Energy  balance  equations 

The  different  elements  of  the  model  are  represented  by  an 
energy  balance  equation,  this  way  the  temperature  field  is 
computed  on  the  whole  surface.  For  the  solids  (i.e.  the  cell 
and  the  flanges)  equations  are  based  on  the  heat  conduction 
for  a  2D  plate  with  source  terms: 

i“"(%i+%i)+i:e  =  0  (3) 

where  the  volumic  heat  sources  gcell  for  the  cell  are  the 
convective  terms,  the  heat  of  the  reaction  and  the  electrical 
power  removed.  For  the  flanges,  the  heat  sources  ^  gfiange 
are  the  radiative  exchange  with  the  oven  and  the  convective 
exchange  with  the  fluids.  The  emissivity  £flange  of  the  flange 


is  estimated  to  be  in  the  range  of  0.4-0. 9  depending  on 
surface  state  [13]. 

For  the  gases,  the  equations  are  similar  for  the  fuel  and  air 
case.  Part  of  the  heat  is  stored  and  removed  by  the  gases  in  the 
form  of  enthalpy  increase.  The  heat  transfer  coefficients  h air 
and  /zfUei  are  computed  from  the  Nusselt  numbers  considering 
the  porous  medium  of  the  zirconia  felt:  this  porous  flow 
enhances  the  heat  transfer.  The  heat  transfer  regime  on  the 
fuel  side  is  dominated  by  the  fuel  conduction  (hydrogen  is 
much  more  conductive  than  the  zirconia  felt),  the  Nusselt 
number  is  then  in  the  order  of  12  [14,15].  On  the  air  side,  the 
Biot  number  has  to  be  computed  first  to  give  a  range  for  the 
Nusselt  number  from  40  to  220.  The  heat  transfer  coefficients 
on  the  air  and  fuel  side  are  consequently  estimated  to  be  2800- 
15,500  and  5750  W/m2  K,  respectively.  Due  to  the  uncer¬ 
tainty  of  the  first  parameter,  a  sensitivity  analysis  is  useful. 

The  equations  for  the  gases  correspond  to  the  local  energy 
balance  of  the  gas  (with  a  sum  on  the  species  i  of  the 
mixture).  The  equation  is: 


where  Y7,  ffpt,  vx ,  L  are  the  molar  concentration,  the  molar 
heat  capacity,  the  gas-velocity  in  the  direction  v  computed 
with  the  potential  flow  model  and  the  thickness  of  the 
zirconia  felt.  A  boundary  condition  is  imposed  at  the  center 
of  the  cell  (gas  inlet)  and  is  computed  with  a  simple  model 
describing  the  conduction  in  the  gas  feed  tube  and  pre¬ 
heating  of  the  gas  before  reaching  the  cell.  The  flow-rates  are 
quite  low  and  in  a  few  centimeters,  the  gases  can  heat-up 
from  oven  temperature  to  the  gas  entry  temperature  (which 
is  close  to  the  cell  temperature  at  the  inlet). 


3.3.  Sensitivity  analysis 


Before  optimizing  the  thermal  parameters,  a  sensitivity 
analysis  is  necessary.  First,  if  the  sensitivity  of  a  parameter  is 
very  low,  it  will  be  difficult  for  the  algorithm  to  find  an 
optimum.  Second,  the  size  of  the  problem  is  reduced.  In 
order  to  select  the  parameters  of  the  thermal  model  to  be 
estimated,  the  sensitivity  analysis  was  done  and  the  results 
summarized  in  Table  1. 

From  the  table,  it  seems  that  the  thermal  conductivity  of 
the  cell  and  the  radiative  emissivity  are  the  most  sensible 
parameters  (considering  the  response  of  the  active  surface 
temperature  Tact  and  the  current). 

All  simulations  were  done  under  the  same  conditions 
except  the  changing  parameter.  They  show  clearly  that 
optimization  has  to  be  performed  on  the  parameters  £flange 
and  2Ceii.  Optimization  on  the  other  parameters  may  lead  to 
problem  because  of  absence  of  sensitivity.  On  Fig.  3,  the 
temperature  field  of  a  cell  operating  at  1.5  A/cm2  is  shown. 
The  temperature  is  maximum  on  the  center  were  reactions 
take  place  and  drops  with  a  very  steep  slope.  On  the  edges, 
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Table  1 


Sensitivity  analysis  for  some  thermal  parameters,  the  values  gives  the 
variation  of  the  results  for  each  parameter 


Parameter 

Sensitivity 

Range 

Tact  (K) 

Current  (A) 

^flange 

0.4-1 

7 

0.1 

Z^fuel 

1000-8000  W/(m2  K) 

0 

0 

h 

“air 

2500-15000  W/(m2  K) 

0.5 

0.01 

■/-cell 

2-9 

2.5 

0.05 

980  mV  at  700  °C.  Nevertheless,  power  density  of  1  W/cm2 
was  achieved  at  700  °C  with  current  density  above  2  A/cm2. 

4.2.  Options  for  an  electrochemical  model 


The  first  trials  used  the  simplest  model  considering  only 
an  ohmic  resistance  depending  upon  temperature.  It  became 
rapidly  obvious  that  this  option  lead  to  large  errors,  espe¬ 
cially  at  low  current.  At  low  current,  activation  losses  cannot 
be  neglected  and  a  new  formulation  was  tested,  taking  into 
account  the  ohmic  losses  and  activation  overpotential.  The 
general  electrochemical  model  gives: 


f^Nernst 


-AG°  RZ  ( (poffm \ 
2  F  2  F  (  pHlo  ) 


U cell  —  f^Nernst  ^  ohmic/'  7  act 


where,  U^emst  is  the  Nemst  potential  and  rjact  the  activation 
overpotential.  As  concentration  polarization  was  never 
observed  even  at  high  current  density,  that  term  was 
neglected.  The  formulations  tested  for  the  activation  over¬ 
potential  were  the  following: 


hsLCt  ho 


2RTsom  .  u 

''”  =  ^u^smh 


-1 


Fig.  3.  Temperature  profile  for  the  cell  at  1.5  A/cm2. 

temperatures  are  quite  high  but  due  to  the  post-combustion 
and  the  low  heat  conductivity,  this  is  of  little  importance  for 
the  temperature  on  the  active  surface. 

4.  Results  and  model  validation 

This  section  presents  the  cell  tests,  results  and  calibration 
of  the  model. 


Eq.  (8),  equivalent  to  Butler- Volmer  behavior,  is  applied  to 
both  anode  and  cathode.  The  exchange  current  j0  is  a 
function  of  temperature.  The  more  empirical  Eq.  (7)  gives 
the  activation  overpotential  as  a  function  of  ri0  and  yp,  which 
are,  respectively,  the  maximum  overpotential  and  the  current 
density  at  which  the  overpotential  reaches  a  plateau.  Ohmic 
losses,  7/act  and  jv  are  functions  of  temperature.  Both  for¬ 
mulations  were  tested  and  parameter  estimation  was  run 
with  the  experimental  data. 

4.3.  Solution  procedure 


4.1.  Cell  tests  and  results 

The  first  set  of  experimental  data  is  from  an  anode 
supported  electrolyte  (ASE)  cell  with  an  LaSrMnCE-YSZ 
(LSM)  composite  cathode  [16]  and  the  second  set  with  an 
ASE  cell  taken  from  a  recent  batch  with  a  LaSrMn03  (LSC) 
cathode  (without  cathode  firing).  The  first  set  of  data  gave  a 
current  density  up  to  1.2  A/cm2  whereas  our  recent  state  of 
the  art  cells  [7]  deliver  up  to  2  A/cm2  current  with  a  LSM- 
LSC  cathode.  The  new  experiment  was  therefore  necessary 
to  measure  higher  current. 

The  conditions  of  the  test  were  flows  of  150  and  250  ml/ 
min  for  hydrogen  and  air,  respectively,  the  oven  was  set  at 
constant  temperature  during  each  current-potential  mea¬ 
surement.  As  LSC  cathode  degradation  is  fast,  measure¬ 
ments  were  done  only  once.  Due  to  pin-holes  on  the 
electrolyte  and  open  circuit  voltage  (OCV)  was  only  around 


In  a  first  step,  optimization  for  each  temperature  is  carried 
out  with  constant  //0,  jp,  A^hmic  and  the  thermal  parameters 
(cell  conductivity  2ceii  and  flange  emissivity  6flange).  The 
gathered  results  were  then  used  to  give  an  initial  guess  of  the 
function  of  temperature  for  each  kinetic  parameter  (which  is 
a  second  order  polynomial  function).  Finally,  optimization  is 
performed  with  the  complete  set  of  data  to  adjust  the 
functions.  In  the  end,  a  continuous  function  for  the  kinetic 
parameters  with  temperature  is  obtained.  The  thermal  para¬ 
meters  seem  to  depend  little  on  temperature  [12,13].  This 
allows  to  keep  them  constant  on  the  whole  temperature 
range  Fig.  4. 

4.4.  Results 

The  model  using  the  simple  Butler- Volmer  equation 
could  not  fit  adequately  the  experimental  data.  Experimental 
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data  may  be  insufficient  to  determine  these  parameters,  or 
the  model  is  not  complete  and  therefore  is  unable  to  describe 
the  physics.  The  second  model  (with  Eq.  (7))  gives  satis¬ 
factory  fits,  Figs.  4  and  5  show  the  current  overpoten¬ 
tial  characteristics  reproduced.  The  current-temperature 


response  (see  Figs.  6  and  7)  is  also  reproduced  with  reason¬ 
able  accuracy,  though  some  under  evaluation  of  the  tem¬ 
perature  remains  at  high  current  for  the  FSM  cathode  case. 

The  offset  remaining  between  the  experiment  and 
simulated  data  (see  on  Fig.  8)  on  the  current-potential 


Fig.  5.  Current-overpotential  characteristics  simulated  and  measured  on  an  ASE  with  LSC  cathode. 
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characteristics  is  caused  by  the  difference  between  the 
experimental  OCV  and  the  simulated  OCV,  which  is  on  the 
order  of  120  mV.  As  the  data  from  the  current  potential 
measurement  was  used  as  an  overpotential  versus  current 
characteristic  for  the  parameter  identification,  the  OCV 
was  not  taken  in  account.  For  the  temperature  response 


(see  Figs.  6  and  7),  the  differences  between  simulation  and 
experiments  are  in  the  order  of  7  K  at  1.2  A/cm2  for  the 
two  high  temperature  measurements  in  the  LSM  case.  One 
explanation  could  be  the  difference  of  OCV,  that  makes 
the  simulation  over-estimate  the  electric  power  removed 
from  the  system. 


Fig.  7.  Temperatures  simulated  and  measured  on  an  ASE  with  LSC  cathode. 
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Fig.  8.  Current-potential  relations  simulated  and  measured  on  an  ASE  with  LSM  cathode. 


Table  2 


Thermal  parameters 

found  for  the  two  cases  by  the  optimization 

/Iceii  (W/(m  K)) 

^flange 

LSM  case 

7.4 

0.36 

LSC  case 

10.1 

0.9 

4.5.  Further  work 

As  the  two  experiments  were  not  performed  in  the  same 
set-up,  the  thermal  parameters  found  for  both  cases  are  quite 
different  (see  Table  2).  To  make  the  routine  electrochemical 
experiment  useful,  the  experiment  with  two  thermocouples 
in  the  same  set-up  and  the  parameter  identification  has  to  be 
repeated.  This  will  ensure  that  the  thermal  parameters  2cen 
and  £fiange  are  constant  for  this  set-up.  Once  enough  con¬ 
fidence  on  those  parameters  is  obtained,  the  model  can  be 
use  to  find  the  kinetic  parameters  from  the  routine  current- 
potential  measurement. 

Regarding  the  kinetic  model,  improvements  are  needed 
first  on  the  OCV  prediction,  as  simulation  over-estimates 
OCV.  An  effort  should  be  made  in  the  direction  of  using  an 
electrochemical  model  to  describe  the  reaction. 

5.  Conclusion 

A  thermal  and  kinetic  model,  representing  a  small  solid 
oxide  fuel  cell  in  a  set-up  used  for  electrochemical  tests  has 
been  developed.  Using  a  parameter  estimation  tool,  the 
simulations  have  been  compared  to  the  experimental  results 


and  gave  satisfactory  fits  for  both  the  current-overpotential 
and  temperature-current  characteristics.  The  remaining 
error  is  mainly  due  to  the  open  circuit  voltage  difference 
between  simulation  and  experiment.  The  measurable  local 
increase  of  temperature  on  the  active  area  can  reach  25-30  K 
at  high  current  depending  on  the  conditions,  as  predicted  and 
confirmed  by  the  model.  Further  experiments  should  be 
carried  out  to  allow  using  the  routine  electrochemical  test 
for  kinetic  parameters  identification.  The  kinetic  model  is 
now  being  implemented  in  a  repeat  element  and  stack  model 
which  is  currently  under  development. 
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